function [elements,processInfo] = InitializeSolutionStep...
    (elements,model,processInfo)
ne = model.countElements();
dim = 2;
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    nn = f.countModes();
    B = zeros(3,dim*nn);
    ng = e.getNumGaussPoint();
    for j=1:ng
        xi = e.getGaussPoint(j-1);
        B = e.computeB(xi(1),xi(2),f,B);
        epsilon = B*ue;
        if strcmp(processInfo.integrationType,'implicit')
            %compute stress
            C = elements(i).constitutiveMatrix;
            elements(i).stress(:,j) = C*epsilon;
        elseif strcmp(processInfo.integrationType,'implex')
            %perform first explicit update
            %compute extrapolated strain-like variable
            kappan = elements(i).kappa_old(j);
            kappanm1 = elements(i).kappa_old_old(j);
            if processInfo.Dt_old == 0
                tmp = 0;
            else
                tmp = processInfo.Dt/processInfo.Dt_old;
            end
            if tmp < 0
                tmp = 0;
            end
            kappa = kappan + tmp*(kappan - kappanm1);
            %compute damage variable based on extrapolated strain-like
            %variable
            d = computeDamageVariable(elements(i),...
                kappa);
            %compute extrapolated stress
            Ce = elements(i).elasticConstitutiveMatrix;
            sigmae = Ce*epsilon;
            elements(i).stress(:,j) = (1-d)*sigmae;
            %compute extrapolated algorithmic tangent
            elements(i).constitutiveMatrix(:,:,j) =...
                        (1-d)*Ce;
        end
    end
end
end
